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abstract 

A general one-dimensional model for the steady adiabatic motion of liquid-volatile mixtures in vertical ducts with 
varying cross-section is presented. The liquid contains a dissolved part of the volatile and is assumed to be incompressible 
and in thermomechanical equilibrium with a perfect gas phase, which is generated by the exsolution of the same volatile. 
An inverse problem approach is used - the pressure along the duct is set as an input datum, and the other physical 
quantities are obtained as output. This fluid-dynamic model is intended as an approximate description of magma- volatile 
mixture flows of interest to geophysics and planetary sciences. It is implemented as a symbolic code, where each line 
stands for an analytic expression, whether algebraic or differential, which is managed by the software kernel independently 
of the numerical value of each variable. The code is versatile and user-friendly and permits to check the consequences of 
different hypotheses even through its early steps. Only the last step of the code leads to an ODE problem, which is then 
solved by standard methods. In the present work, the model is applied to the study of two sample cases, representing 
the ascent of magma-gas mixtures on the Earth and on a Jupiters satellite, Io. Both cases lead to approximate but 
realistic descriptions of explosive eruptions, by taking pressure curves as inputs and outputting conduit shapes together 
£jrj with mixture density, temperature and velocity along the ducts. 
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1 Introduction 

A number of authors have modeled conduit flows of magma- volatile mixtures in order to represent the physical 
phenomena involved in volcanic eruptions. These flows are studied in the cases of steady effusive and explosive 
volcanic activity (MC GETCHIN k ULLRICH 1973, WILSON et alii 1980, WILSON & HEAD 1981, KIEFFER 
1982, SLEZIN 2003), by using either of two conduit models to set boundary conditions in order that the 
system might be solved semi-analytically. More specifically, conduits are taken as being either parallel-sided 
or lithostatically pressure-balanced, the former model being characterized by walls parallel to one another and 
perfectly rigid (WILSON et alii 1980, WILSON & HEAD 1981), the latter to a scenario where local changes in 
flow pressure relative to lithostatic pressure in the country rocks are promptly accommodated by wall failure 
(VALENTINE & GROVES, 1996). A consequence of the parallel-sided scenario is the occurrence of the transonic 
transition of flow regime (M = 1) at the surface, whereas the lithostatically pressure-balanced mode implies 
wallrock deformation and erosion, and assumes that the conduit is capable of adjusting itself such that stress 
across the walls is zero. In the pressure-balanced scenario, the flow regime transition is allowed to migrate 
downwards, from the surface to some depth in the upper conduit. As a result, upper-conduit velocities can 
be supersonic, but shocks and rarefaction waves may form at or after the transition zone (KIEFFER 1982, 
CHAPMAN 2000), which may cause massive changes in velocity, temperature and pressure in the fluid passing 
through them. Importantly, shock generation is very likely to be associated with non-gradual changes in conduit 
profile and/or abrupt variations in flow pressure gradients. 

The present work is concerned with steady pressure-balanced scenarios, where the conduit shape evolution 
has come to an end and the conduit profile doesn't change with time. The presented model deals mainly with 
explosive volcanic activity, but allows for effusive activity as a subcase. We model the ascent of magma and 
magma-gas mixtures from the base of the conduit up to the vent region. The magma rising through the crust 
is allowed to incorporate variable amounts of crustal volatile deposits. The volatile can be partly dissolved and 
partly exsolved (gas phase) in the magma, and for gas proportions larger than 10 wt%, the mixture behaves 
as an ideal pseudogas (WALLIS 1969, KIEFFER 1982, LU & KIEFFER 2009). Other model assumptions 
are taken from modelers of adiabatic flow of magma-gas mixtures through volcanic conduits (KIEFFER 1982, 
BURESTI & CASAROSA 1989 and 1993, MASTIN 1995 and 2000), where no heat is transferred across the 
conduit walls during the eruption. The supersonic velocities resulting from the gain in kinetic energy occurring 
in pressure-balanced scenarios are counter-balanced by losses in elastic and thermodynamic energy. 
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In the present model, the flow is one dimensional, steady and homogenous so, at any given depth, flow 
properties can be averaged across the entire cross-sectional area of the conduit. The gas phase, which is taken 
to be incondensable, behaves essentially as an ideal gas whereas the liquid phase is incompressible and in 
thermodynamic and mechanical equilibrium with the gas. Finally, at each depth the pressure of the mixture is 
taken to be uniform and equal to that of the gas phase alone. The model is completely analytic until it leads to 
an ordinary differential equation (or a set of ODEs in the most general case), which is then solved numerically by 
standard methods. This is obtained by suitable hypotheses, which simplify the theoretical treatment maintaining 
the capability of yielding predictions for the main physical quantities involved. These hypotheses, exposed in 
detail in the next sections, can be summarized as follows: the liquid phase (magma+dissolved volatile) is 
Newtonian; in this phase the dissolved mass fraction of volatile is in general much lower than the magma mass 
fraction; the exsolution process as a function of depth is slow; the conduit shape and the mixture temperature 
vary gradually; the mixture velocity exhibits large variations only in the upper, final part of the conduit. 

The corresponding author has written a Mathematica© symbolic code - its advantage lying in its versatile 
tailorability for specific research needs and goals. Thanks to the symbolic approach, each line of the code is 
analytic like the corresponding equation of model, up to the last step, where a discrete computational component 
is still present, in that numerical methods for the solution of an ODE are used in computations. In this way, 
checks about the modification of an hypothesis, which may potentially give rise to a different model, can be 
easily done due to the fact that, as the code is running, many substitutions of an equation into a part of another, 
as well as derivations, integrations and so forth are being performed by the software. However, it is worth to 
remember that the appearance of non-physical or inconsistent calculations shall always be prevented by the 
author of the code. 

Among the possible applications, two case studies of activity on the Earth and Jupiters satellite, Io, are 
henceforth considered, though the analysis could be equally well-suited for any other planetary body in the 
Solar System for which evidence of present and/or past volcanic activity exists. 

2 Fundamental equations 

It is considered an unidimensional flow of a magma+volatile mixture in a vertical conduit, where the liquid 
fraction (magma+dissolved volatile) is Newtonian and the exsolved volatile is a perfect gas, according to the 
hypotheses presented in the introduction above. The conduit is extended over the vertical coordinate domain 
zo < z < 0, where z = is the conduit upper end at surface level (vent) and zq is the depth of the conduit 
lower end, the location of the latter being referred to as a virtual reservoir. This means that it may represent 
a real magma+volatile reservoir, or a zone where magma coming from deeper zones enters the conduit, with a 
chance of interacting with a given amount of crustal volatile deposits. The conduit cross section A is considered 
circular in what follows, A(z) — itr 2 {z) where r(z) is the conduit radius, but this hypothesis is not so restrictive 
and can be easily generalized to different shapes. Similarly, the generalization from vertical to oblique conduits 
is not difficult, and it can be obtained by projecting the gravity acceleration on the direction of interest. It 
is assumed that the variation of the conduit cross section is gradual, i.e. the radius varies slowly. A formal 
criterion for this property will be introduced in ij3.4l and 13.51 

The equation set to consider in this problem is made of mass, momentum and energy conservation laws for 
a magma-gas mixture, and is written as follows: 



where the prime denotes z-derivatives. Here v(z),p(z), p(z),T(z) are velocity, pressure, density and temperature 
of the mixture, g is the gravitational acceleration and / is a friction factor which can be related to the wall 
properties of the conduit (SCHLICHTING 1979). / is typically constant for turbulent flows, but it could become 
z-dependent for low Reynolds numbers. The energy equation involves also the constant pressure specific heat 
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of the mixture c p (z), the magma mass fraction in the mixture m(z) and the magma density p m . In the most 
general case, all quantities depend on z except g and p m . The state equation of the mixture is written as 



1 -mp/p m 1 

where R = n(z)R g is a generalized form of the gas constant R g , accounting for the gas mass fraction n(z) — 
1 — m(z). It represents the pressure equilibrium between magma and gas, and the denominator accounts for the 
volumetric fraction of magma. For large gas mass fractions, the state equation behaves as the pseudogas state 
formula p = pRT+ small corrections, which in turn tends to the ideal gas state law for n —> 1. Formally, this 
state equation and relations ([TJ . . . ([3]) form a system of 4 equations that characterize the problem under study. 

Different approaches to the problem are possible. Considering <?, /, p m , m(z), c p (z) as parameters of the 
model, to assign directly or to relate to the unknown functions by means of constitutive relationships, there are 
5 unknown functions, namely v(z),p(z), p{z), T(z) and r(z). Starting from a formal system of 4 equations, two 
methods can be considered: 

• Direct problem solution: the conduit shape is known by hypothesis together with the accessory conditions 
(function values at reference points, like the reservoir depth). The method leads to the flow pressure along the 
conduit and to the other functions v{z), p(z), T(z), .... 

• Inverse problem solution: the flow pressure in the conduit p(z) is known by hypothesis together with the 
accessory conditions. The method leads to the conduit shape as a radius curve r(z) and to the other functions 
characterizing the flow. 

The present work deals with an inverse problem, simplified by suitable assumptions, and implemented as a 
symbolic code, as disclosed in the introduction. A remarkable part of the analytic work is performed when the 
code is running until the desired result is achieved, i.e. until a final equation or set of equations is obtained, 
which is then solved numerically in the same code by means of high-level instructions. 



3 Model operation 

3.1 Starting data and lower conduit laws 

All formulas in what follows correspond to lines of the corresponding symbolic code. At first, as an input, 
a pressure curve p(z) must be reasonably defined in a physically meaningful way on the vertical coordinate 
domain zq < z < 0. Then, since the physics is remarkably influenced by the gaseous phase (exsolved volatile) 
present in the mixture, it must be considered that the present model deals with two main kinds of flow: if at 
a given depth Zf the gas volume fraction Vf of exsolved volatile contained in the mixture overcomes a critical 
value which is generally believed to be close to 0.75 (point of fragmentation), then the mixture velocity will 
increase dramatically, thus entering a compressible flow regime similar to that of ideal gases, with a chance for 
supersonic flow conditions to develop in the upper conduit and a potential to generate an explosive eruption. 
Alternatively, an effusive event will occur. 

In the available literature, more than one critical value of Vf is reported, which could all well suit the present 
model; in contrast, the formulation proposed here is not suitable for implementing fragmentation criteria based 
on the extensional-strain rate of a mixture element dv/dz. Most importantly, in each particular case under 
study, it is not even known in advance whether or not fragmentation will occur at a certain z = Zf, which may 
be considered as a boundary region separating conduit portions characterised by different flow regimes. For this 
reason, the vertical domain must be considered in general as the union of two distinct domains and the pressure 
formally defined as the union of two separate curves, 

^ IPa(z) for zo < z < Zf (lower conduit) 
lp{,(z) for Zf < z < (upper conduit) . 

Here p(z) is considered continuous across the whole domain zq < z < 0, but a discontinuity in its first derivative 
at Zf is allowed. This property may also appear in other relationships of this model without spoiling the validity 
of results, since the final equations can be solved over the two domains with suitable matching conditions. On 
the other hand, the definition (0 may also be written in the form p(z) — 6{zf — z)p a (z) + 6{z — Zf)pt(z) 
where 0(z) is the step function, and the model can be easily generalized introducing a smoothing taper function 
s(z— Zf), such as an hyperbolic tangent based on the conduit radius scale, and working with a pressure definition 
of the kind p{z) = s{zf — z)p a (z) + s(z — Zf)pb(z). In the simplest cases, fragmentation doesn't take place and 
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the "lower conduit" coincides with the overall domain. Examples of initial guesses for pressure in realistic cases 
will be given in 

In this model, the fragmentation depth and pressure are obtained to a certain degree of accuracy. The total 
volatile mass fraction nt must be set as an input datum, and may involve a number of volatile species, as follows: 



E 



nu- (6) 



Each volatile species can be dissolved/exsolved in/from the magma in different amounts. For each species, 
models of dissolved fractions are found in the literature, and they all lead to the total fraction of dissolved 
volatile, as follows: 

n d {z) = Y t n di {z) = Y t F t {z,...) (7) 

i i 

where the symbols Fi represent different dissolution models, which typically depend on p(z) but may also affect 
one another (examples will be given in Sj4]). Once the total dissolved fraction is known, the exsolved volatile in 
the magma, i.e. the gas mass fraction of the mixture, can be written as 

n(z) =n t -n d (z), (8) 

and then the liquid (magma+dissolved volatile) fraction is m(z) = 1 — n(z). The density pi of this liquid should 
account for the dissolved volatile mass fraction in the liquid, n v , and the dissolved volatile density, p v , but the 
relevant formula can be simplified as follows, under the general assumptions n v -C 1 and p v <C p m '- 

Pi = T\ \ i - P™- ( 9 ) 

(1 - n v )p v + n v p m 

In this model the variations of n ( i(z),n(z) and m(z) are taken to be very small and gradual, in order to 
neglect the relevant derivatives, a condition satisfied in many real cases. In a perturbative formulation based on 
a small parameter e, this property would be expressed for example as n = n(ez) = n(Z) so that dn/dz = e dn/dZ 
and the calculations would be performed at order 0(1) by neglecting terms of order 0(e). 

Also, the total specific heat of any volatile species can be obtained by combining the heats of the different 
species (again, examples will be given in Q: 

CraO) = -^YlM^Cpgi, ( 10 ) 

where rii(z) — nu — ndi(z) is the exsolved fraction for the i-th species, and the specific heats c pg i are taken as 
mean values over a short range of temperatures around the lower end value To. This approximation relies on 
the hypothesis that, over a very wide range of eruption conditions studied in the literature, the temperature 
differences between the initial and final sections of the ducts are generally found to be limited to a few tens of 
degrees; in other words, the present model refers to flows where AT /To <C 1. In an analogous way, the specific 
gas constant for any gaseous volatile can be obtained as: 

n \ z ) i 

where the R g i are the gas constants for the species involved. Now, using the relation c vg — c pg — R g and 
accounting for the specific heat of the magma c m (input datum), the specific heats of the mixture at constant 
pressure and volume can be written as 

c p (z) = n(z)c pg (z) + m(z)c m (12) 
c v (z) = n(z)c vg (z) + m(z)c m . (13) 

As a consequence of the previous hypotheses, also c p ,c v and R are slowly varying parameters. 

To determine the fragmentation, it is necessary to know the gas density p g , which is given here by the ideal 
law written for the lower conduit 

Pga = Tr|r; (14) 
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this formula involves the temperature T a in the lower conduit, which is expressed in this model by integrating 
over z the energy conservation ([3]), starting from the lower end values Tq,pq,vq at z = zq (input data), and 

j 2 ) are small 



simplifying the result under the hypothesis that in the lower conduit the velocity variations (v, 
in comparison with the other terms in the original equation: 



T a (z) = T + 



1 



1 



Pn 



(PO - Pa) + 77O0 - V ) + 9 ( z - z ) 



T - 
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(po - Pa) +g(z Q - z) 



(15) 



Of course, the physical meaning of this approximation is that in the lower conduit the radius changes very 
slowly and the gas volume is small, which is generally true for a wide class of flows of this kind; the effects of 
this hypothesis can be checked as explained in §3.51 Then, the gas volume fraction can be written as 



Vf = 



np„ 



np r , 



mpga 



(16) 



where the contribution of the dissolved volatile was neglected after equation ([9]). Finally, assuming that frag- 
mentation takes place for a known critical value Vf c (typically 0.75), and solving the equation Vf(zf) = Vf c , 
the fragmentation depth Zf is obtained. The relevant pressure is pf — p a {zf) = p (zf), so that the two domains 
(upper and lower conduits) are completely determined. Examples of initial guesses for pressure, as well as a 
check of the effects of the approximations introduced in this section will be given in 



3.2 Temperature and density on the complete domain 



Up to now, the temperature and density of the mixture have only been defined in the lower conduit. The 
treatment of the upper conduit in the most general case could turn out to be much more complicated, but 
in the present work it is simplified by focusing on the main physics involved in this region. An approximated 
and useful model for the temperature can be defined by considering that, after the fragmentation point Zf, 
the dominant phenomenon is mixture expansion, and more specifically the expansion of the gas fraction as the 
magma is disrupted into small particles. Under the hypothesis of thermodynamic equilibrium, an isentropic 
expansion law for the mixture can be obtained as in BUREST1 & CASAROSA 1989, combining the isentropic 
condition ds — c p dT/T — R dp/p — with the state law of the mixture ((4]), and considering the z-derivatives of 
c p and R negligible as in the previous section. The resulting law can be expressed as 



or 



' 1 - m p/p m \ 7 
p = const., 



1-7 

p i T = const., 



where 7 is the ratio of the specific heats of mixture at constant pressure and volume (|T2l and (fT3")l , 

7(2) = c p (z)/c v (z). 

The considerations above suggest expressing the temperature in the following way: 
\T a (z) = T + 1 

T = 




-ipo - Pa) + g(z Q - Z) 



Ta(Zf) 



Ph_ 

Pf 



for Zq < Z < Zf 

for Zf < z < 0, 



(17) 
(18) 
(19) 

(20) 



where the lower conduit law follows relation (|15p and the upper conduit law follows relation (|18l) representing 
an expansion starting from the value T a {zf) and determined by the pressure variation Pb{z) up to the final or 
vent pressure value pf (input data). Even if simplified, expression (|2"0|) for temperature doesn't prevent the 
appearance of phenomena like the adiabatic heating of the mixture, as will be shown in S|4j Like p{z), the 
definition (|2U)) may have a discontinuity in its first derivative, which doesn't affect the problem solution. Also 
in this case, it can be removed by tapering. 

The gas density p g in the upper conduit is still given by the ideal law, here p g b = Pb/ (R g Tb), so that the 
global density can be written as: 



PgPn 



mpg + np m 



(21) 
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As the gas volume fraction (I16p . this formula is simplified by neglecting the very small contribution of the 
dissolved volatile. Also in p(z), a first derivative discontinuity at z = zt may be present, because it involves p 
and T. This discontinuity is removable by tapering, and the same considerations as above hold for this function 
and all of the followings, when they involve p and/or p, T. 

3.3 Velocity and Mach number, complete domain 

The flow speed can be related to the other functions by means of the continuity equation ([TJ . The velocity law 
turns out to be 

p{z)r 2 (z) 

where the input radius, velocity and density 7"o,«o and po are known as input data at the conduit lower end. 
At this stage v(z) is formally unknown, because r(z) is unknown. This is not a problem for the symbolic code, 
which retains the analytic expression of v in the computer's memory. 

The sound speed formula for the mixture, thanks to the hypothesis of thermodynamic equilibrium, can be 
obtained by relation (fTT|) using the definition c 2 = (dp/dp) s and the state equation Q: 

(23) 

1 - mp/p m 

where j(z) is given by (|19|) and R — n(z)R g . As in the state equation, the denominator accounts for the 



volumetric fraction of magma. For large gas mass fractions, formula (|23|) expresses the sound speed in a 
pseudogas and finally, for n 1, in an ideal gas. 

The Mach number can be defined as usual, i.e. M(z) = v(z)/c(z). As the velocity, M(z) is at this stage 
known only in symbolic form. It is worth noting that, if supersonic flow develops, the sonic coordinate along 
the conduit can no longer be determined by the ideal condition dr/dz — 0, as shown in the next subsection. 

3.4 Final integration 

Up to now, continuity and energy conservation as well as relations having the physical meaning of a state law 
have been used in the model build-up. On the other hand, the momentum equation ||2} relates all the quantities 
defined up to now and it can be used for the problem solution. If the variations of p and p are related to the 
speed of sound through relation p' = (dp/dp) s p' = p' /c 2 , the momentum equation can be rewritten as 

Since the present method is completely analytic up to the integration of equation (|24[) , by substitution of the 
previous formulas in this equation, a huge expression is obtained, which is far too cumbersome for reading. This 
expression is stored in the computer's memory by the symbolic code, and after substitution of the scalar input 
parameters, it can be solved for r(z) by standard methods on the two domains Zo < z < Zf and Zf < z < with 
the boundary condition r(z ) = r and the continuity condition at z = Zf (the continuity condition is no longer 
necessary if the pressure definition (|S|) has a continuous derivative over the whole domain). In the symbolic 
code all steps (lines) but the last one are analytic. The last one instead provides the numerical solution, which 
is obtained through a single, high level, final instruction. The relevant numerical method can be automatically 
selected by the Mathematica© software among the standard ones (Runge-Kutta, Adams...) in order to reach 
the highest level of accuracy in compliance with the internal representation of numbers in the computer used, 
and the integration step is selected in an adaptive way with the same goal. Alternatively, an experienced user 
can drive the solver in a more detailed way. 

Equation (pM)) shows also that when M = 1 the conduit radius satisfies the condition 



r 4+f 



where the right hand side is always positive, so that the sonic condition can be reached only at positions where 
r' > 0. This equation collapses into the ideal one if g and / vanish. 

The radius obtained by integrating equation (I24[) must vary gradually in order to satisfy the basic hypotheses, 
and this can be expressed by the condition that the radius variation over a length scale in the order of the same 
radius must be small: 

\r(z + Az) — r(z)\ < r(z) for Az~r(z) and z < z < -r(0). (25) 



G 



3.5 Accuracy checks 



The effect of the approximations introduced in the model can be checked a posteriori in different ways. A 
possible method relies on the use of energy conservation in integral form, 

e t = c p T + — p + \-v 2 + gz (26) 

Pm ^ 

which relates the total energy per unit mass e t to pressure, velocity and temperature curves. In particular, at 
the lower end depth, equation (|2"B1) reads 

e = c p0 To + Po + §«g + gz . 
After the solution of the problem, i.e. the integration of equation (|24|) . the relevant energy et can be calculated 
by substituting in equation (1261) the quantities c p ,T, m,p and v determined by the present model, and the 
resulting function e*(z) can be compared with the reference value eo- For a "perfect model" et should remain 
constant along z, i.e. et = bq. In the present case the ratio 

fo^z) (2?) 
eo 

becomes an index of the method accuracy, representing the relative error on the total energy. As a general 
criterion, results could be validated if this index remains less than a given threshold, well below the unity. In 
fact, this quantity appears to be particularly sensitive to non-physical temperature variations: if the present 
model and code are used within the range set by the original hypotheses, the value of (|27|) amounts to less 
than a few percentage points, otherwise it grows rapidly. Similar indices can be defined for other quantities of 
physical interest. 

A different but important validity check can be obtained by putting condition (1251) in a linearized differential 
form, which gives rise to the quantity 

\r'(z)\, (28) 

i.e. a non-dimensional expression for radius variation rapidity. In order to have a gradually varying duct, |r'(z)| 
should remain well below unity along the conduit, and a criterion of comparison with a given threshold can 
be easily introduced. Another significant quantity is \r' [z)/r{z)\, which has the dimensions of a wavenumber. 
It may represents the scale of conduit shape fluctuations along z and it should remain well below 1/r. Large 
values of |r'(z)| distributed over long z ranges may indicate a violation of unidimensionality, whereas large local 
variations of \r'(z)/r(z)\ may represent a violation of the iscntropic assumption, such as that implied by the 
appearance of shocks in the compressible part of the flow. 



4 Sample results 

4.1 Simplified formulas for flow pressure 

For the sake of testing the method presented in section S}3j some examples of "pressure hypothesis" are given 
here. In the simplest case, making the basic assumptions of incompressible flow in the subsonic range, the 
pressure could be taken as a linearly decreasing function; under conditions of constant mass fractions and heat 
coefficients, and slowly varying gas volume fraction, the same assumption would lead to slowly varying radii 
and velocities along the conduit, close to the extreme of pipe flow with friction where in the standard model r 
and v are constants. Thus, the simplest nontrivial hypothesis for the pressure is 

/ s Pe-PQ, v 

P{z) =Po (z ~ Zq) 

where the input data are the magma driving pressure po of steady eruptions of magma and magma-gas mixtures 
(for which a plausible range of values can be chosen from available models) and the output pressure p e (known 
by models and/or direct observations). 

However, the increase of the gas volume fraction has important effects - if and when at a given depth zt the 
volume of exsolved volatile in the mixture overcomes the critical value and magma fragmentation occurs, even 
a pressure gradient which is linear at large depths could give rise to a compressible flow on a duct. In general, 
when a flow of this kind enters the compressible regime, there is also a deviation of the overall pressure curve 
from the linear behaviour; this deviation is often not negligible, even if it remains small in many interesting 
cases. Soon after magma fragmentation, the pressure typically decreases following a curve which is remarkably 



7 



zo z f 

Figure 1: Sketch of real pressures compared with approximated pressures. Solid lines, thick and thin: examples of real 
pressure. Dashed lines: piecewise linearized pressures. The z extent of the nonlinear zone in this figure is shown enlarged 
for clarity. 
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Figure 2: Particular of the approximated pressure curve, close-up in the fragmentation zone. Thick line: example of 
real pressure as in fig. [1] Dashed line: piecewise linearized pressure. 



nonlinear only in the very upper conduit, this trend being qualitatively shown in figure [TJ where the two solid 
curves (thick and thin) show different kinds of nonlinear behaviours. The z extent of the nonlinear zone in this 
figure must be shown in an enlarged view for clarity, since often is very small with respect to the overall length 
of the conduit. For a wide class of eruption flows, in the lower conduit, which represents the most vertically 
extensive portion of the duct, and according to WILSON & HEAD (1981), a real pressure curve in the lower 
conduit is expected to be similar to the lithostatic pressure in the surrounding rocks because otherwise the 
conduit would be destroyed or heavily modified. The lithostatic pressure is typically very close to the linear 
behaviour or weakly nonlinear. The property of Zf to be very close to the end of the conduit, i.e. \zf \ <C \zq\, 
implies that the local deviation of p(z) from a simple linear behaviour is not large, and in general the real flow 
pressure, a piecewise linearized approximation and the real lithostatic pressure are all close together, with the 
maximum differences becoming apparent only over a short (upper) portion of the conduit. 

For these reasons, here it is considered a piecewise linearization of p(z) which introduces only a small error 
on the overall behaviour. The figure [1] shows that in the following examples the real pressure p(z), which can 
be considered as made of two parts pa(z) for z < zt and pb{z) for z > Zf (solid lines), can be approximated 
by piecewise linearized pressures p a (z) and Pb(z) in the same subdomains (dashed lines). Here the input data 
are still Po,p e whereas the fragmentation depth and pressure Zf and pf can be calculated as in £13.11 with small 
errors. The details of real and approximated pressure curves for the thick line of figure [T] in the fragmentation 
zone are sketched in figure [5J 

Formally, the approximated pressure curve proposed for the following examples is a piecewise linear function 
on the vertical coordinate domain zq < z < 0: 

{Pa(z) = po H —(z — zq) for zq < z < Zf (lower conduit) 

Pe-Pf ' . ( 29 ) 

Pb(z) = Pf -{z — Zf) for Zf < z < (upper conduit) 
z f 

The 1st derivative discontinuity can be avoided by tapering, as described in M3 - 1 1 In equation (|29[) , po is the 
driving pressure at the lower end (depth Zq), which can be taken to be similar to the value of the lithostatic 
pressure for a fixed depth, with possible corrections accounting for varying crustal and tectonic settings across 
the planets/satellites under scrutiny; p e is the output pressure at z = 0, whilst Ap is an estimation of the 
pressure drop along a conduit having a length in the order of L = \z \. This estimation can be obtained from 



S 



the known law for pipes with friction, namely 

&P = -f P U 2 ^ 5 (30) 

where / is the friction factor, p a scale density, which can be taken at this step equal to the magma density 
p m , U a scale velocity, in the order of the initial velocity parameter vq, D a scale diameter which can be taken 
equal to 2ro (twice the initial radius of the conduit). 

A code starting from equation (|29j) can be implemented by defining as a first step p a (z) and following the 
procedure of M3.1l to determine the fragmentation depth Zf. Then pb(z) can be defined, and all the steps of H3.2I 
to 13.41 can be followed. If zj turns out to be or greater, fragmentation doesn't take place, the lower conduit 
coincides with the whole duct, and equation (|29|) collapses into the linear expression for p a (z). This is a model 
that suits an effusive eruption, in this case the input datum p e (exit or vent pressure) is not used, and the output 
pressure is automatically determined as the value of p Q (0). This doesn't prevent the appearance of variations 
of radius and velocity, owing to the variation of gas volume fraction, but in general they are remarkably smaller 
than the variations in a flow with fragmentation. 

Finally, it is worth to remember that the present model for pressure is just one of the possible and viable 
choices, and the symbolic code permits the exploration of different hypotheses in an easy fashion. 

4.2 An example based on terrestrial data. 

As a first example, it is considered a mixture of terrestrial magma with CO2 and H2O. To set the volatile 
solubility properties, the magma is assumed to be basaltic, since basalt is the most common rock type found 
on the terrestrial planets. The input data at the conduit lower end are reported in what follows: 

z = -5000m ; T = 1350K ; p = 7 ■ 10 7 Pa ; p m = 2900kg/m 3 ; r = 2m ; v = 25 m/s . 
In this example the input velocity represents a fully developed flow more than a flow starting from a real reservoir. 
The total mass fractions of the volatile species are ntco 2 =0.01 ; n t H 2 o =0.02 . The volatile properties in this 
kind of magma are modeled after GERLACH (1986) and MASTIN (1995): the exsolved volatile mass fraction 
nt is given by equation l[5|). which reads here 

n(z) = n C02 {z) + hh 2 o(z) = [™tco 2 ~ n dC o 2 ( z )} + [nm 2 ~ n d n 2 o(z)}- 
The dissolved fractions for the two species are calculated in different ways, and in particular ndco 2 I s a linear 
function of the input pressure 

ndco 2 ( z ) = 01 + a-2P{z), 

whereas ridH 2 o depends on ndco 2 and must be obtained considering an algebraic system of the kind 

n d n 2 o = F(p,p p ) 
ndH 2 o = «tH 2 o - G{n dC o 2 ,p,p p ) 

where the complicated functions F and G are reported in the cited works (GERLACH 1986, MASTIN 1995) and 
also depend on the partial pressure p p of H2O in the mixture. The solution of the system gives nda 2 o(z) and 
p P (z). Then, the heat coefficients are obtained from (|10| and (fTTjl . The remaining input data are g = — 9.81m/s 2 ; 
p e = 10 5 Pa; / = 0.06. The value of friction factor / represents a high roughness; here it is important to 
recall that the present definition of / is consistent with the pressure drop expressed by equation (|30p . as in 
SCHLICHTING (1979), but other authors use definitions which differ by a factor of four from the present one. 

The input pressure curve modeled as in H4.1\ is shown in figure [3] and the results are shown in figures 2] to 
[7] A tapering was applied only to the input function p(z). Output (vent) values of the physical quantities are 
reported in the captions for this particular case. The global error on the energy conservation after equation (|27p 
is 2.4%. The non-dimensional steepness after equation (|28p is less than 0.024 up to 5 radii below the vent, and 
reaches the value 0.125 in the last radius below the vent. The scenario depicted in this example could represent 
an active lava fountain or the Plinian phase of an explosive eruption on the Earth. 

It can be seen that the modification of few initial data, like the pressure profile, the total volatile mass fraction 
or the output pressure, leads to remarkable variations of the vent data such output velocity and divergence angle, 
which could represent different terrestrial eruptions. 

4.3 An example based on data for satellite Io. 

Another example, chosen to represent some properties of the eruptions on Jupiters satellite, Io, is presented in 
what follows. In the relevant literature, the magmatic materials are believed to be basaltic and/or ultramafics 
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Figure 3: Terrestrial example, input pressure curve: overall a and close-up in the upper zone b. A weak knee is visible at 
the fragmentation depth, which is in this case -209m ~ 0.04[«o |. The 1st derivative discontinuity is removed by tapering 
over a 50 diameters scale. 
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Figure 4: Terrestrial example, density a and temperature b. The density exhibits a progressive decrease as the gas 
volume fraction increases. As said in the model description, the temperature is slightly raised by subsonic adiabatic 
heating, then decreases in the final expansion. The overall temperature variation ratio is AT/Tq ~ 18%. 





V 




[ra/s 


z [ra] 



-5000 -4000 -3000 -2000 -1000 



100 
80 
60 
40 
20 



a) 



z [m] 



5000 



4000 



3000 



2000 



1000 



b) 



Figure 5: Terrestrial example, velocity a and Mach number b. The flow is sonic (M = 1) at z = —34.4m. Output 
velocity is 95.8m/s~ 4«o and output Mach number is 1.81. This value of M is related to the low value of the sound 
speed, which is in the mixture remarkably lower than in surrounding air. 
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Figure 6: Terrestrial example, conduit diameter (not to scale). White level is effectively proportional to gas volume 
fraction. 
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Figure 7: Terrestrial example, conduit diameter, close-up in the output zone, scaled coordinates. White level propor- 
tional to gas volume fraction. The output radius is 6.97m ~ 3.5ro and the final semidivergence angle is 18.0 degrees. 
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Figure 8: Example on Io, input pressure curve: overall a and close-up in the upper zone b. The fragmentation depth 
is in this case -774m ~ 0.026|zq|- The 1st derivative discontinuity is removed by tapering over a 75 diameters scale. 



in composition, and the most important volatile species is assumed to be SO2. This information is used here 
to set such necessary parameters, as the magma density and the volatile solubility properties. The input data 
at the conduit lower end are listed below: 

z = -30000m ; T = 1700K ; p = 1.2 • 10 8 Pa ; p m = 3000kg/m 3 ; r = 0.42m ; v Q = 6 m/s . 
In this example the input velocity may represent the starting speed of the mixture after the initial development 
length at the reservoir-conduit connection. The total volatile mass fraction is n t — 0.015, which is likely to 
partly arise from the incorporation of gaseous/liquid SO2 from shallow crustal aquifers, as shown in previous 
modelling of ionian eruptions (CATALDO et alii 2002). These crustal deposits are taken to be in proximity to 
zq. The volatile properties are modeled after MYSEN 1977, expressing the dissolved fraction rid(z) as a power 
series on p of the kind 

n d{z) = J2j a jP 3 ^ 

approximated by a 3-d order truncation. The remaining input data are g = — 1.8m/s 2 ; p e = 2 • 10 4 Pa; / = 0.055. 

The input pressure curve is shown in figure [5] and the results are shown in figures [5] to 1121 A tapering was 
applied only to the input function p(z). The global error on the energy conservation after equation (|27p is in 
this case 1.9%. The non-dimensional steepness after equation (f28|) is less than 0.027 up to 5 radii below the 
vent, and becomes 0.099 at the last radius. 




Figure 9: Example on Io, density a and temperature b. The overall temperature variation ratio is AT /To ~ 10%. 
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Figure 10: Example on Io, velocity a and Mach number b in the final part of the conduit. The flow is sonic at 
z — —11.9m. Output velocity is 89.9m/s~ 14vo and output Mach number is 1.57. This value, which is related to the 
exit pressure chosen as an input datum, suggests the appearance of an underexpanded jet at the conduit vent since the 
surrounding ambient has a very low pressure. This is likely to represent the lowermost portion of a small-scale erupting 
plume on Io. 
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Figure 11: Example on Io, conduit diameter (not to scale). White level is effectively proportional to gas volume fraction. 
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Figure 12: Example on Io, conduit diameter, close-up in the output zone, scaled coordinates. White level proportional 
to gas volume fraction. The output radius is 2.37m ~ 6ro and the final semidivergence angle is 10.5 degrees. 
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5 Conclusions 



The present paper deals with the study of the motion of magma- volatile mixtures in long vertical ducts with 
friction, originated from the desire to develop a model for the flow of magmatic fluid along volcanic conduits 
during steady phases of explosive eruptions, and allowing for effusive eruptions as a subcase. A simplified 
approach to this problem was taken, i.e. the mixture was assumed to be composed of an incompressible liquid 
phase and a perfect gas phase in conditions of thermomechanical equilibrium. The flow was treated as one- 
dimensional, homogeneous, steady and adiabatic. By means of further simplifying hypotheses, a "quasi-analytic" 
model capable of yielding predictions for the main physical quantities involved was obtained. In particular, the 
model is completely analytic through the whole series of steps except for the final ordinary differential equation 
(which becomes a set of ODEs in the most general case), that is solved numerically by standard methods. 

The model is implemented as a symbolic code, where each line is similar to an equation of the model itself, or 
to what a human would write on a blackboard. Formal substitution of one equation into another and many other 
steps like derivations and integrations are performed by the software at running time, including the solution of 
the final ODE. The code is easy to modify, offering the choice to change at any time the basic hypotheses of 
the model and/or introducing new ideas. The degree of fulfilment of the assumptions upon which the model is 
based is checked by calculating suitable indicators a posteriori. 

The model permits rapid analyses of the effects of the variation of the input parameters, a brief example 
was included for the sake of clarity, but without the aim of giving an exhaustive treatment, that could be the 
object of a further work. 

Many further improvements in the model may be envisaged, like the study of non-circular cross sections, 
the introduction of different state equations, that could for example represent volatiles present in the mixture 
as supercritical fluids, and also the modification of the conservation laws in order to represent the interaction 
of a developed mixture flow with a crustal deposit of volatile that adds mass to the mixture in a given range of 
depths. Again, these improvements could be considered as objects of future works. 
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